; docformat = 'rst'
;
;+
;
; :Purpose:
;   Filter data using Kolmogorov-Zurbenko filter
;   
; :Inputs:
;   x_in: array of locations of the data to be filtered (e.g. time)
;   y_in: array to be filtered
;   sigma: array of uncertainties in y_in
;   
; :Keywords:
;   iterations: (input, integer) numeber of filter passes
;   window: (input, floating) width of filter window in the same units as x_in
;   growth: (input, Boolean) set to true to filter the growth rate of y_in
;   min_elements: (input, integer), minimum number of data points that must lie within a particular 
;     window in order for average to be calculated
;   out_window: (output, floating) effective 'window' size of smoothed data
;   
; :Outputs:
;   Structure containing:
;     x: Array of locations of smoothed data
;     y: Array of smoothed data, or growth rate (if /growth set)
;
; :Example::
;   IDL> x=findgen(100)*0.1
;   IDL> y=sin(x) + randomn(10, 100)
;   IDL> out=mr_kz_filter(x, y, window=0.5, iterations=4, out_window=out_window)
;   IDL> print, out_window
;         1.00000
;   IDL> plot, x, y
;   IDL> oplot, out.x, out,y
;   
; :History:
; 	Written by: Matt Rigby, University of Bristol, Nov 8, 2012
;
;-
function mr_kz_filter, x_in, y_in, growth=growth, iterations=iterations, $
  window=window, min_elements=min_elements, out_window=out_window

  on_error, 2

  if keyword_set(min_elements) eq 0 then begin
    min_elements=1
  endif

  if keyword_set(growth) then begin
    y=(y_in[1:-1] - y_in[0:-2])/(x_in[1:-1] - x_in[0:-2])
    x=(x_in[0:-2] + x_in[1:-1])/2.
  endif else begin
    y=y_in
    x=x_in
  endelse
  
  ys=y
  xStart=x[0] + window
  xEnd=x[-1] - window
  xiStart=(where(x ge xStart))[0]
  xiEnd=(where(x gt xEnd))[0]
  
  xs=x
  
  for k=0, iterations-1 do begin
    ys_prev=ys
    xs_prev=xs
    for xi=0, n_elements(xs)-1 do begin
      wh=where(xs_prev gt xs_prev[xi]-window/2. and $
        xs_prev le xs_prev[xi]+window/2. and $
        finite(ys_prev), count)
      if count ge 1 then begin
        ys[xi]=mean(ys_prev[wh], /nan)
      endif else begin
        ys[xi]=!values.d_nan
      endelse
    endfor
  endfor
  
  xs[0:xiStart]=!values.f_nan
  xs[xiEnd:-1]=!values.f_nan
  
  ys[0:xiStart]=!values.f_nan
  ys[xiEnd:-1]=!values.f_nan
  
  out_window=window*sqrt(float(iterations))
  
  return, {x:x, y:ys}

end